Exploration bulkRNAseq analyses

Accompanying ‘CONVOCALS_downstream’

General Setup

Create a new analysis directories.

- general directory

- stains directory

- for plots

- for output of summary results

- for baseline tables

source("scripts/functions.R")

source("scripts/pack02.packages.R")

* General packages...

* Genomic packages...

Today = format(as.Date(as.POSIXlt(Sys.time())), "%Y%m%d")
Today.Report = format(as.Date(as.POSIXlt(Sys.time())), "%A, %B %d, %Y")
source("scripts/colors.R")

This notebook

In this notebook we explore the files for the bulk RNAseq analyses and select the subset of data we need.

Background

Here we obtain data from the CONVOCALS_downstream in plaques.

• Genes.xlsx - list of genes of interest.
• Variants.xlsx - list of variant(s) of interest.

library(openxlsx)

gene_list_df <- read.xlsx(paste0(TARGET_loc, "/targets.xlsx"), sheet = "Genes")

gene_list <- unlist(gene_list_df$Gene)
gene_list

 [1] "SCGB3A2"  "LIX1"     "IGSF9B"   "ND6"      "ALB"      "OR10A4"   "RASEF"    "NEDD4"    "TRIM49D1" "TCL1A"    "ND4L"     "ATP8"     "FBXO15"   "F5"       "TMEM212" 
[16] "PTPRD"    "CYP46A1"  "OR2T33"   "SORBS2"   "ITGA7"    "RNASE1"   "FOS"      "HMOX1"    "LAPTM5"   "MMP9"     "SIGLEC1"  "FTL"      "CD14"     "HCST"     "TIMP3"   
[31] "CCL2"     "SAT1"     "CD163"    "PTGDS"    "LGALS9"   "ACKR1"    "NT5DC2"   "TGFBI"    "C1QC"     "S100A9"  

variant_list <- read.xlsx(paste0(TARGET_loc, "/targets.xlsx"), sheet = "Variants")

DT::datatable(variant_list)

NA

Load data

First we will load the data:

• prepared bulk RNA sequencing (RNAseq) experimental data from carotid plaques
• Athero-Express clinical data.

Now you need to choose which of the datasets you will use for the downstream analyses.

• raw
  • not transformed, raw counts
• normalized data
  • normalized based on size factors based on estimateSizeFactors()
• VST data
  • normalized based on size factors and log-transformed using vst()

Note: we advise using vst() transformed normalized data

Log-transformed normalized data

AERNAScomboSE <- readRDS(file = paste0(OUT_loc, "/",Today,".AERNAScomboSEvst.CEA.1092pts.SE.after_qc.IC_academic.RDS"))
# AERNAScomboSE <- readRDS(file = paste0(OUT_loc, "/20240109.AERNAScomboSEvst.CEA.1092pts.SE.after_qc.IC_academic.RDS"))

Clinical Data

First, we extract the clinical data.

bulkRNA_meta_clin <- as.tibble(colData(AERNAScomboSE))
bulkRNA_meta_clin %<>%
  # mutate(macrophages = factor(macrophages, levels = c("no staining", "minor staining", "moderate staining", "heavy staining"))) %>%
  # mutate(smc = factor(smc, levels = c("no staining", "minor staining", "moderate staining", "heavy staining"))) %>%
  # mutate(calcification = factor(calcification, levels = c("no staining", "minor staining", "moderate staining", "heavy staining"))) %>%
  # mutate(collagen = factor(collagen, levels = c("no staining", "minor staining", "moderate staining", "heavy staining"))) %>%
  # mutate(fat = factor(fat, levels = c("no fat", "< 40% fat", "> 40% fat"))) %>%
  mutate(study_number_row = STUDY_NUMBER) %>%
  as.data.frame() %>%
  column_to_rownames("study_number_row")
names(bulkRNA_meta_clin)[names(bulkRNA_meta_clin) == "STUDY_NUMBER"] <- "study_number"
head(bulkRNA_meta_clin)

dim(bulkRNA_meta_clin)

[1] 1092   69

Second, we define the variables we need.

# Baseline table variables
basetable_vars = c("Hospital", "ORyear", "Artery_summary",
                   "Age", "Gender",
                   # "TC_finalCU", "LDL_finalCU", "HDL_finalCU", "TG_finalCU",
                   "TC_final", "LDL_final", "HDL_final", "TG_final",
                   # "hsCRP_plasma",
                   "systolic", "diastoli", "GFR_MDRD", "BMI",
                   "KDOQI", "BMI_WHO",
                   "SmokerStatus", "AlcoholUse",
                   "DiabetesStatus",
                   "Hypertension.selfreport", "Hypertension.selfreportdrug", "Hypertension.composite", "Hypertension.drugs",
                   "Med.anticoagulants", "Med.all.antiplatelet", "Med.Statin.LLD",
                   "Stroke_Dx", "sympt", "Symptoms.5G", "AsymptSympt", "AsymptSympt2G",
                   "Symptoms.Update2G", "Symptoms.Update3G",
                   "restenos", "stenose", "Artery_summary",
                   "CAD_history", "PAOD", "Peripheral.interv",
                   "EP_composite", "EP_composite_time", "epcom.3years", 
                   "EP_major", "EP_major_time","epmajor.3years",
                   "MAC_rankNorm", "SMC_rankNorm", "Macrophages.bin", "SMC.bin",
                   "Neutrophils_rankNorm", "MastCells_rankNorm",
                   "IPH.bin", "VesselDensity_rankNorm",
                   "Calc.bin", "Collagen.bin",
                   "Fat.bin_10", "Fat.bin_40",
                   "OverallPlaquePhenotype", "Plaque_Vulnerability_Index")

basetable_bin = c("Gender", "Artery_summary",
                  "KDOQI", "BMI_WHO",
                  "SmokerStatus", "AlcoholUse",
                  "DiabetesStatus",
                  "Hypertension.selfreport", "Hypertension.selfreportdrug", "Hypertension.composite", "Hypertension.drugs",
                  "Med.anticoagulants", "Med.all.antiplatelet", "Med.Statin.LLD",
                  "Stroke_Dx", "sympt", "Symptoms.5G", "AsymptSympt", "AsymptSympt2G",
                  "Symptoms.Update2G", "Symptoms.Update3G",
                  "restenos", "stenose", "Artery_summary",
                  "CAD_history", "PAOD", "Peripheral.interv",
                  "EP_composite", "Macrophages.bin", "SMC.bin",
                  "IPH.bin",
                  "Calc.bin", "Collagen.bin",
                  "Fat.bin_10", "Fat.bin_40",
                  "OverallPlaquePhenotype", "Plaque_Vulnerability_Index")
# basetable_bin

basetable_con = basetable_vars[!basetable_vars %in% basetable_bin]
# basetable_con

Expression differences

From here we can analyze whether specific genes differ between groups, or do this for the entire gene set as part of DE analysis, and then select our genes of interest.

From there, we extract the gene expression values, plus gene identifier, and select the genes of our interest.

expression_data <- assay(AERNAScomboSE)

# extract expression values from vsd, including ensembl names
expression_data <- as_tibble(data.frame(gene_ensembl = rowRanges(AERNAScomboSE)$feature_id, assay(AERNAScomboSE))) %>%
     mutate_at(vars(c("gene_ensembl")), list(as.character)) ## gene_ensembl needs to be character for annotation to work

mutate_at: no changes

# annotations
# gene symbol - via org.Hs.eg.db
# columns(org.Hs.eg.db)
expression_data$symbol <- mapIds(org.Hs.eg.db,
                    keys = expression_data$gene_ensembl,
                    column = "SYMBOL",
                    keytype = "ENSEMBL",
                    multiVals = "first")

'select()' returned 1:many mapping between keys and columns

# tidy and subset
expression_data_sel <- expression_data %>%
     dplyr::select(gene_ensembl, symbol, everything()) %>%
     # filter(symbol == "APOE" | symbol == "TRIB3") %>% # filter APOE and TRIB3
     dplyr::filter(symbol %in% gene_list)

head(expression_data_sel)

# tidy and subset non-selected genes
set.seed(141619)
expression_data_sample <- expression_data %>%
     dplyr::select(gene_ensembl, symbol, everything()) %>%
     sample_n(1000) %>%
     unite(symbol_ensembl, symbol, gene_ensembl, sep = "_", remove = FALSE)

sample_n: removed 20,835 rows (95%), 1,000 rows remaining

expression_data_sample_mean <- expression_data_sample %>%
  select_if(is.numeric) %>%
  colMeans() %>%
  as_tibble(rownames = "study_number") %>%
  dplyr::rename(expression_value_sample = value)

select_if: dropped 3 variables (symbol_ensembl, gene_ensembl, symbol)

Furthermore, the expression_data_sel df was gathered into a long form df for annotation with symptoms variables from the vsd object, and later visualization and statistics.

# gather expression_data_sel df into long df form for annotation, plotting and statistics
expression_long <-
     gather(expression_data_sel, key = "study_number", value = "expression_value", -c(gene_ensembl, symbol))

gather: reorganized (ae1, ae1026, ae1029, ae1032, ae1054, …) into (study_number, expression_value) [was 40x1094, now 43680x4]

# old school way
# Annotate with smoking variables
# sample_ids <- expression_long$study_number
# mm <- match(expression_long$study_number, rownames(colData(vsd)))
#
# ## Add traits to df
# ## Binary traits
# expression_long$sex <- colData(vsd)$sex[mm]
# expression_long$testosterone <- colData(vsd)$testosterone[mm]
# expression_long$t_e2_ratio <- colData(vsd)$t_e2_ratio[mm]

# new school way
plaque_phenotypes_cat <- c("Macrophages.bin",
                           "SMC.bin",
                           "Calc.bin",
                           "Collagen.bin",
                           "Fat.bin_10", 
                           # "Fat.bin_40",
                           "IPH.bin")

plaque_phenotypes_num <- c("MAC_rankNorm", #"macmean0",
                           "SMC_rankNorm", #"smcmean0",
                           # "MastCells_rankNorm", #"mast_cells_plaque",
                           # "Neutrophils_rankNorm", #"neutrophils",
                           "VesselDensity_rankNorm")

expression_long <- expression_long %>%
  left_join(bulkRNA_meta_clin %>% dplyr::select(study_number,
                                       plaque_phenotypes_cat,
                                       plaque_phenotypes_num,
                                       epmajor.3years, 
                                       #epmajor.30days,
                                       AsymptSympt2G,
                                       Gender, Hospital, StudyName),
            by = "study_number") %>%
  mutate(epmajor_3years_yn = str_replace_all(epmajor.3years, c("Excluded" = "yes", "Included" = "no"))) #%>%

left_join: added 14 columns (Macrophages.bin, SMC.bin, Calc.bin, Collagen.bin, Fat.bin_10, …)           > rows only in x                               0           > rows only in bulkRNA_meta_clin %>% d.. (     0)           > matched rows                            43,680           >                                        ========           > rows total                              43,680

  # mutate(epmajor.30days_yn = str_replace_all(epmajor.30days, c("Excluded" = "yes", "Included" = "no")))

head(expression_long)

# expression_long %>%
#   write_tsv("genes_interest_expression.txt")

Gene expression - distribution

Filter genes

In case some genes are not available in our data we could filter them here.

gene_list

 [1] "SCGB3A2"  "LIX1"     "IGSF9B"   "ND6"      "ALB"      "OR10A4"   "RASEF"    "NEDD4"    "TRIM49D1" "TCL1A"    "ND4L"     "ATP8"     "FBXO15"   "F5"       "TMEM212" 
[16] "PTPRD"    "CYP46A1"  "OR2T33"   "SORBS2"   "ITGA7"    "RNASE1"   "FOS"      "HMOX1"    "LAPTM5"   "MMP9"     "SIGLEC1"  "FTL"      "CD14"     "HCST"     "TIMP3"   
[31] "CCL2"     "SAT1"     "CD163"    "PTGDS"    "LGALS9"   "ACKR1"    "NT5DC2"   "TGFBI"    "C1QC"     "S100A9"  

This code is just an example to filter the list from genes that are not in the data.

• COL3A ==> not found
• COL2A ==> not found

gene_list_rm <- c("COL3A", "COL2A") 

temp = gene_list[!gene_list %in% gene_list_rm]

gene_list_qc <- c(temp)

# gene_list_qc <- gene_list
# 
# for debug
# gene_list_qc_replace <- c("MRTFA")

Plotting expression

Figure 1: Expression of genes of interest: boxplots

# Make directory for plots
ifelse(!dir.exists(file.path(QC_loc, "/Boxplots")), 
       dir.create(file.path(QC_loc, "/Boxplots")), 
       FALSE)

[1] FALSE

BOX_loc = paste0(QC_loc,"/Boxplots")

for(GENE in gene_list_qc){
  cat(paste0("Plotting expression for ", GENE,".\n"))
  temp <- subset(expression_long, symbol == GENE)
  
  compare_means(expression_value ~ Gender, data = temp)
  p1 <- ggpubr::ggboxplot(temp,
                          x = "Gender",
                          y = "expression_value",
                          color = "Gender",
                          palette = "npg",
                          add = "jitter",
                          ylab = paste0("normalized expression ", GENE,"" ),
                          repel = TRUE
                          ) + stat_compare_means()
  print(p1)
  cat(paste0("Saving image for ", GENE,".\n"))
  
  ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_gender.png"), plot = last_plot())
  ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_gender.pdf"), plot = last_plot())
  ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_gender.ps"), plot = last_plot())
  
  compare_means(expression_value ~ Gender, data = temp)
  p2 <- ggpubr::ggboxplot(temp,
                          x = "Gender",
                          y = "expression_value",
                          color = "StudyName",
                          palette = "npg",
                          add = "jitter",
                          ylab = paste0("normalized expression ", GENE,"" ),
                          repel = TRUE
                          ) #+ stat_compare_means()
  print(p2)
  cat(paste0("Saving image for ", GENE,".\n"))
  
  ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_Gender_by_AERNAStudy.png"), plot = last_plot())
  ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_Gender_by_AERNAStudy.pdf"), plot = last_plot())
  ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_Gender_by_AERNAStudy.ps"), plot = last_plot())
  
  compare_means(expression_value ~ StudyName, data = temp)
  p3 <- ggpubr::ggboxplot(temp,
                          x = "StudyName",
                          y = "expression_value",
                          color = "StudyName",
                          palette = "npg",
                          add = "jitter",
                          ylab = paste0("normalized expression ", GENE,"" ),
                          repel = TRUE
                          ) + stat_compare_means()
  print(p3)
  cat(paste0("Saving image for ", GENE,".\n"))
  
  ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_AERNAStudy.png"), plot = last_plot())
  ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_AERNAStudy.pdf"), plot = last_plot())
  ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_AERNAStudy.ps"), plot = last_plot())
  
  compare_means(expression_value ~ StudyName, data = temp)
  p4 <- ggpubr::ggboxplot(temp,
                          x = "StudyName",
                          y = "expression_value",
                          color = "Gender",
                          palette = "npg",
                          add = "jitter",
                          ylab = paste0("normalized expression ", GENE,"" ),
                          repel = TRUE
                          ) #+ stat_compare_means()
  print(p4)
  cat(paste0("Saving image for ", GENE,".\n"))
  
  ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_AERNAStudy_by_Gender.png"), plot = last_plot())
  ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_AERNAStudy_by_Gender.pdf"), plot = last_plot())
  ggsave(filename = paste0(BOX_loc, "/", Today, ".",GENE,".expression_vs_AERNAStudy_by_Gender.ps"), plot = last_plot())
  rm(temp, p1, p2, p3, p4 )
}

Plotting expression for SCGB3A2.
Saving image for SCGB3A2.
Saving image for SCGB3A2.
Saving image for SCGB3A2.
Saving image for SCGB3A2.
Plotting expression for LIX1.
Saving image for LIX1.
Saving image for LIX1.
Saving image for LIX1.
Saving image for LIX1.
Plotting expression for IGSF9B.
Saving image for IGSF9B.
Saving image for IGSF9B.
Saving image for IGSF9B.
Saving image for IGSF9B.
Plotting expression for ND6.
Saving image for ND6.
Saving image for ND6.
Saving image for ND6.
Saving image for ND6.
Plotting expression for ALB.
Saving image for ALB.
Saving image for ALB.
Saving image for ALB.
Saving image for ALB.
Plotting expression for OR10A4.
Saving image for OR10A4.
Saving image for OR10A4.
Saving image for OR10A4.
Saving image for OR10A4.
Plotting expression for RASEF.
Saving image for RASEF.
Saving image for RASEF.
Saving image for RASEF.
Saving image for RASEF.
Plotting expression for NEDD4.
Saving image for NEDD4.
Saving image for NEDD4.
Saving image for NEDD4.
Saving image for NEDD4.
Plotting expression for TRIM49D1.
Saving image for TRIM49D1.
Saving image for TRIM49D1.
Saving image for TRIM49D1.
Saving image for TRIM49D1.
Plotting expression for TCL1A.
Saving image for TCL1A.
Saving image for TCL1A.
Saving image for TCL1A.
Saving image for TCL1A.
Plotting expression for ND4L.
Saving image for ND4L.
Saving image for ND4L.
Saving image for ND4L.
Saving image for ND4L.
Plotting expression for ATP8.
Saving image for ATP8.
Saving image for ATP8.
Saving image for ATP8.
Saving image for ATP8.
Plotting expression for FBXO15.
Saving image for FBXO15.
Saving image for FBXO15.
Saving image for FBXO15.
Saving image for FBXO15.
Plotting expression for F5.
Saving image for F5.
Saving image for F5.
Saving image for F5.
Saving image for F5.
Plotting expression for TMEM212.
Saving image for TMEM212.
Saving image for TMEM212.
Saving image for TMEM212.
Saving image for TMEM212.
Plotting expression for PTPRD.
Saving image for PTPRD.
Saving image for PTPRD.
Saving image for PTPRD.
Saving image for PTPRD.
Plotting expression for CYP46A1.
Saving image for CYP46A1.
Saving image for CYP46A1.
Saving image for CYP46A1.
Saving image for CYP46A1.
Plotting expression for OR2T33.
Saving image for OR2T33.
Saving image for OR2T33.
Saving image for OR2T33.
Saving image for OR2T33.
Plotting expression for SORBS2.
Saving image for SORBS2.
Saving image for SORBS2.
Saving image for SORBS2.
Saving image for SORBS2.
Plotting expression for ITGA7.
Saving image for ITGA7.
Saving image for ITGA7.
Saving image for ITGA7.
Saving image for ITGA7.
Plotting expression for RNASE1.
Saving image for RNASE1.
Saving image for RNASE1.
Saving image for RNASE1.
Saving image for RNASE1.
Plotting expression for FOS.
Saving image for FOS.
Saving image for FOS.
Saving image for FOS.
Saving image for FOS.
Plotting expression for HMOX1.
Saving image for HMOX1.
Saving image for HMOX1.
Saving image for HMOX1.
Saving image for HMOX1.
Plotting expression for LAPTM5.
Saving image for LAPTM5.
Saving image for LAPTM5.
Saving image for LAPTM5.
Saving image for LAPTM5.
Plotting expression for MMP9.
Saving image for MMP9.
Saving image for MMP9.
Saving image for MMP9.
Saving image for MMP9.
Plotting expression for SIGLEC1.
Saving image for SIGLEC1.
Saving image for SIGLEC1.
Saving image for SIGLEC1.
Saving image for SIGLEC1.
Plotting expression for FTL.
Saving image for FTL.
Saving image for FTL.
Saving image for FTL.
Saving image for FTL.
Plotting expression for CD14.
Saving image for CD14.
Saving image for CD14.
Saving image for CD14.
Saving image for CD14.
Plotting expression for HCST.
Saving image for HCST.
Saving image for HCST.
Saving image for HCST.
Saving image for HCST.
Plotting expression for TIMP3.
Saving image for TIMP3.
Saving image for TIMP3.
Saving image for TIMP3.
Saving image for TIMP3.
Plotting expression for CCL2.
Saving image for CCL2.
Saving image for CCL2.
Saving image for CCL2.
Saving image for CCL2.
Plotting expression for SAT1.
Saving image for SAT1.
Saving image for SAT1.
Saving image for SAT1.
Saving image for SAT1.
Plotting expression for CD163.
Saving image for CD163.
Saving image for CD163.
Saving image for CD163.
Saving image for CD163.
Plotting expression for PTGDS.
Saving image for PTGDS.
Saving image for PTGDS.
Saving image for PTGDS.
Saving image for PTGDS.
Plotting expression for LGALS9.
Saving image for LGALS9.
Saving image for LGALS9.
Saving image for LGALS9.
Saving image for LGALS9.
Plotting expression for ACKR1.
Saving image for ACKR1.
Saving image for ACKR1.
Saving image for ACKR1.
Saving image for ACKR1.
Plotting expression for NT5DC2.
Saving image for NT5DC2.
Saving image for NT5DC2.
Saving image for NT5DC2.
Saving image for NT5DC2.
Plotting expression for TGFBI.
Saving image for TGFBI.
Saving image for TGFBI.
Saving image for TGFBI.
Saving image for TGFBI.
Plotting expression for C1QC.
Saving image for C1QC.
Saving image for C1QC.
Saving image for C1QC.
Saving image for C1QC.
Plotting expression for S100A9.
Saving image for S100A9.
Saving image for S100A9.
Saving image for S100A9.
Saving image for S100A9.

Figure 2A: Expression of genes of interest: histograms

# Make directory for plots
ifelse(!dir.exists(file.path(QC_loc, "/Histograms")), 
       dir.create(file.path(QC_loc, "/Histograms")), 
       FALSE)

[1] FALSE

HISTOGRAM_loc = paste0(QC_loc,"/Histograms")

for(GENE in gene_list_qc){
  # cat(paste0("Plotting expression for ", GENE,".\n"))
  temp <- subset(expression_long, symbol == GENE)
  p1 <- ggpubr::gghistogram(temp,
                          x = "expression_value",
                          y = "..count..",
                          color = "Gender", fill = "Gender",
                          palette = "npg",
                          add = "median",
                          ylab = paste0("normalized expression ", GENE,"" )  
                          )
  print(p1)
  cat(paste0("Saving image for ", GENE,".\n"))
  ggsave(filename = paste0(HISTOGRAM_loc, "/", Today, ".",GENE,".distribution.png"), plot = last_plot())
  ggsave(filename = paste0(HISTOGRAM_loc, "/", Today, ".",GENE,".distribution.pdf"), plot = last_plot())
  ggsave(filename = paste0(HISTOGRAM_loc, "/", Today, ".",GENE,".distribution.ps"), plot = last_plot())

  rm(temp, p1 )
}

Saving image for SCGB3A2.
Saving image for LIX1.
Saving image for IGSF9B.
Saving image for ND6.
Saving image for ALB.
Saving image for OR10A4.
Saving image for RASEF.
Saving image for NEDD4.
Saving image for TRIM49D1.
Saving image for TCL1A.
Saving image for ND4L.
Saving image for ATP8.
Saving image for FBXO15.
Saving image for F5.
Saving image for TMEM212.
Saving image for PTPRD.
Saving image for CYP46A1.
Saving image for OR2T33.
Saving image for SORBS2.
Saving image for ITGA7.
Saving image for RNASE1.
Saving image for FOS.
Saving image for HMOX1.
Saving image for LAPTM5.
Saving image for MMP9.
Saving image for SIGLEC1.
Saving image for FTL.
Saving image for CD14.
Saving image for HCST.
Saving image for TIMP3.
Saving image for CCL2.
Saving image for SAT1.
Saving image for CD163.
Saving image for PTGDS.
Saving image for LGALS9.
Saving image for ACKR1.
Saving image for NT5DC2.
Saving image for TGFBI.
Saving image for C1QC.
Saving image for S100A9.

Figure 2B: Expression of genes of interest: density plots

# Make directory for plots
ifelse(!dir.exists(file.path(QC_loc, "/Density")), 
       dir.create(file.path(QC_loc, "/Density")), 
       FALSE)

[1] FALSE

DENSITY_loc = paste0(QC_loc,"/Density")

for(GENE in gene_list_qc){
  # cat(paste0("Plotting expression for ", GENE,".\n"))
  temp <- subset(expression_long, symbol == GENE)
  p1 <- ggpubr::gghistogram(temp,
                          x = "expression_value",
                          y = "..density..",
                          color = "Gender", fill = "Gender",
                          palette = "npg",
                          add = "median",
                          ylab = paste0("normalized expression ", GENE,"" )  
                          )
  print(p1)
  cat(paste0("Saving image for ", GENE,".\n"))
  ggsave(filename = paste0(DENSITY_loc, "/", Today, ".",GENE,".density.png"), plot = last_plot())
  ggsave(filename = paste0(DENSITY_loc, "/", Today, ".",GENE,".density.pdf"), plot = last_plot())
  ggsave(filename = paste0(DENSITY_loc, "/", Today, ".",GENE,".density.ps"), plot = last_plot())
  
  rm(temp, p1 )
}

Saving image for SCGB3A2.
Saving image for LIX1.
Saving image for IGSF9B.
Saving image for ND6.
Saving image for ALB.
Saving image for OR10A4.
Saving image for RASEF.
Saving image for NEDD4.
Saving image for TRIM49D1.
Saving image for TCL1A.
Saving image for ND4L.
Saving image for ATP8.
Saving image for FBXO15.
Saving image for F5.
Saving image for TMEM212.
Saving image for PTPRD.
Saving image for CYP46A1.
Saving image for OR2T33.
Saving image for SORBS2.
Saving image for ITGA7.
Saving image for RNASE1.
Saving image for FOS.
Saving image for HMOX1.
Saving image for LAPTM5.
Saving image for MMP9.
Saving image for SIGLEC1.
Saving image for FTL.
Saving image for CD14.
Saving image for HCST.
Saving image for TIMP3.
Saving image for CCL2.
Saving image for SAT1.
Saving image for CD163.
Saving image for PTGDS.
Saving image for LGALS9.
Saving image for ACKR1.
Saving image for NT5DC2.
Saving image for TGFBI.
Saving image for C1QC.
Saving image for S100A9.

Compare expression to the expression of a sample of 1,000 genes

Figure 3: comparing expression of genes of interest to mean expression of a sample of 1,000 random genes

expression_wide <- expression_long %>%
  tidyr::unite("symbol_ensembl", symbol:gene_ensembl, remove = TRUE, sep = "_") %>%
  # dplyr::select(-gene_ensembl, -symbol) %>%
  tidyr::spread(key = symbol_ensembl, value = expression_value) 

# the next 3 lines of code gave an error when selecting for genes_interest, since one of the genes of interest is missing: FGF3 is not in the data set. So, we need to select for the other 15 genes.

temp <- expression_long %>%
  tidyr::unite("symbol_ensembl", symbol:gene_ensembl, remove = TRUE, sep = "_") %>%
  dplyr::select(symbol_ensembl)
temp <- temp %>% 
  tidyr::separate(symbol_ensembl, c("symbol", "gene_ensembl"), remove = FALSE) %>%
  dplyr::distinct(., symbol_ensembl, .keep_all = TRUE)

symbol_ensembl_gene_list_qc <- temp[temp$symbol %in% gene_list_qc,]$symbol_ensembl

expression_wide2 <- expression_wide %>%
  left_join(expression_data_sample_mean, by = "study_number") %>%
  dplyr::select(study_number, StudyName, Gender, symbol_ensembl_gene_list_qc, expression_value_sample)

left_join: added one column (expression_value_sample)           > rows only in x                              0           > rows only in expression_data_sample_.. (    0)           > matched rows                            1,092           >                                        =======           > rows total                              1,092

expression_long2 <- expression_wide2 %>%
  gather(gene, expression_value, -study_number, -StudyName, -Gender) %>%
  mutate(gene = str_replace_all(gene, c("expression_value_sample" = "Random genes"))) %>%
  tidyr::separate(gene, c("gene", "symbol_ensembl"), remove = FALSE, sep = "_") %>%
  dplyr::select(-symbol_ensembl) %>%
  dplyr::mutate(gene = factor(gene, levels = c("Random genes", gene_list_qc)))

gather: reorganized (CYP46A1_ENSG00000036530, NEDD4_ENSG00000069869, IGSF9B_ENSG00000080854, FTL_ENSG00000087086, SIGLEC1_ENSG00000088827, …) into (gene, expression_value) [was 1092x44, now 44772x5]Warning: Expected 2 pieces. Missing pieces filled with `NA` in 1092 rows [43681, 43682, 43683, 43684, 43685, 43686, 43687, 43688, 43689, 43690, 43691, 43692, 43693, 43694, 43695, 43696, 43697, 43698, 43699, 43700, ...].

expression_long2

mean_1000_genes <- mean(expression_data_sample_mean$expression_value_sample)
# head(expression_long2)
#

p1 <- ggpubr::ggboxplot(expression_long2,
                          x = "gene",
                          y = "expression_value",
                          color = uithof_color[16],
                          add = "jitter",
                          add.params = list(size = 3, jitter = 0.3), 
                          ylab = paste0("normalized expression ")
                          ) +
    geom_hline(yintercept = mean_1000_genes, linetype = "dashed", color = uithof_color[26], size = 0.4) + 
    theme(axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1), # change orientation of x-axis labels
          axis.text.y = element_text(size = 12),
          axis.title.x = element_text(size = 14),
          axis.title.y = element_text(size = 14)) 
  p1
  
  ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes.png"), plot = last_plot())

Saving 7.29 x 4.51 in image

  ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes.pdf"), plot = last_plot())

Saving 7.29 x 4.51 in image

  ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes.ps"), plot = last_plot())

Saving 7.29 x 4.51 in image

p2 <- ggpubr::ggboxplot(expression_long2,
                          x = "gene",
                          y = "expression_value",
                          color = "StudyName",
                          palette = "npg",
                          add = "jitter",
                          add.params = list(size = 3, jitter = 0.3), 
                          ylab = paste0("normalized expression ")
                          ) +
    geom_hline(yintercept = mean_1000_genes, linetype = "dashed", color = uithof_color[26], size = 0.4) + 
    theme(axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1), # change orientation of x-axis labels
          axis.text.y = element_text(size = 12),
          axis.title.x = element_text(size = 14),
          axis.title.y = element_text(size = 14)) 
  p2
  
  ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes_by_AERNAStudy.png"), plot = last_plot())

Saving 7.29 x 4.51 in image

  ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes_by_AERNAStudy.pdf"), plot = last_plot())

Saving 7.29 x 4.51 in image

  ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes_by_AERNAStudy.ps"), plot = last_plot())

Saving 7.29 x 4.51 in image

  
  rm(p1, p2, temp )

p3 <- ggpubr::ggboxplot(expression_long2,
                          x = "gene",
                          y = "expression_value",
                          color = "Gender",
                          palette = "npg",
                          add = "jitter",
                          add.params = list(size = 3, jitter = 0.3), 
                          ylab = paste0("normalized expression ")
                          ) +
    geom_hline(yintercept = mean_1000_genes, linetype = "dashed", color = uithof_color[26], size = 0.4) + 
    theme(axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1), # change orientation of x-axis labels
          axis.text.y = element_text(size = 12),
          axis.title.x = element_text(size = 14),
          axis.title.y = element_text(size = 14)) 
  p3
  
  ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes_by_Gender.png"), plot = last_plot())

Saving 7.29 x 4.51 in image

  ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes_by_Gender.pdf"), plot = last_plot())

Saving 7.29 x 4.51 in image

  ggsave(filename = paste0(PLOT_loc, "/", Today, ".TargetExpression_vs_1000genes_by_Gender.ps"), plot = last_plot())

Saving 7.29 x 4.51 in image

  
  rm(p3)

List all the gene for all unique genes where their mean expression is higher than the mean of the 1,000 random genes (mean_1000_genes).

# calculate the mean for each gene in expression_long2
expression_long2_mean <- expression_long2 %>%
  group_by(gene) %>%
  summarise(mean_expression = mean(expression_value))

summarise: now 41 rows and 2 columns, ungrouped

expression_long2_mean

expression_long2_high <- expression_long2_mean %>%
  dplyr::filter(mean_expression > mean_1000_genes) %>%
  dplyr::select(gene) %>%
  dplyr::distinct(gene, .keep_all = TRUE)
expression_long2_high

List all genes with expression lower than mean_1000_genes.

expression_long2_low <- expression_long2_mean %>%
  dplyr::filter(mean_expression < mean_1000_genes) %>%
  dplyr::select(gene) %>%
  dplyr::distinct(gene, .keep_all = TRUE)

expression_long2_low

List all genes where the mean of their expression is higher or lower than the mean of 1,000 random genes, and where there is a statistical difference between gender.

# 1. Calculate mean expression per gender per gene
expression_long2_mean_gender <- expression_long2 %>%
  dplyr::group_by(gene, Gender) %>%
  dplyr::summarise(mean_expression = mean(expression_value), .groups = 'drop')

# Print the mean expression per gender
print(expression_long2_mean_gender)

# 2. Perform a Wilcoxon rank-sum test (Mann-Whitney U test) for gender differences per gene
expression_long2_mean_gender_t_test_results <- expression_long2 %>%
  dplyr::group_by(gene) %>%
  dplyr::summarise(wilcox_test = list(wilcox.test(expression_value ~ Gender)), .groups = 'drop') %>%
  dplyr::mutate(p_value = sapply(wilcox_test, function(x) x$p.value)) %>%
  dplyr::filter(p_value <= 0.05) %>% # Filter for p-values below 0.05
  dplyr::arrange(p_value)  # Order by p-value in ascending order

# Print the p-values from the t-tests
print(expression_long2_mean_gender_t_test_results)

# 3. Write the results to an Excel file using openxlsx
# Create a new workbook
wb <- createWorkbook()

# Add a worksheet
addWorksheet(wb, "Significant Genes")

# Write the dataframe to the worksheet
writeData(wb, "Significant Genes", expression_long2_mean_gender_t_test_results)

# Save the workbook to an Excel file
saveWorkbook(wb, paste0(OUT_loc, "/", Today,".iph_area_075_vst_dgea_gender_significant_genes.xlsx"), overwrite = TRUE)

Heatmaps for genes of interest

If we would put these correlations in one simple and comprehensible figure, we could use a correlation heatmap. Again, correlation coefficients used here are Spearman’s.

Figure 4: correlation heatmap between expression of genes of interest

library(tidyverse)
library(magrittr)

temp <- expression_wide %>%
  column_to_rownames("study_number") %>%
  dplyr::select(symbol_ensembl_gene_list_qc) %>% 
  drop_na() %>% # drop NA 
  Filter(function(x) sd(x) != 0, .) # filter variables with sd = 0

# average out columns that start with the same gene symbol, where the gene symbol is the first part of the column name separated by '_'; rename the columns to the part before '_'
# Load necessary library
library(dplyr)

# Assume your data frame is named `expression_wide`

# Step 1: Identify the columns that correspond to gene symbols using regex
gene_columns <- grep("ENSG\\d+", colnames(temp), value = TRUE)

# Step 2: Extract gene symbols using regular expression
# Extract everything before the underscore (_) as the gene symbol
gene_symbols <- sub("(_ENSG\\d+)", "", gene_columns)

# Step 3: Create a list of unique gene symbols
unique_genes <- unique(gene_symbols)

# Step 4: For each unique gene, either average the columns (if multiple) or rename (if single)
temp2 <- temp

for (gene in unique_genes) {
  # Find all columns that match the gene symbol
  gene_cols <- grep(paste0("^", gene, "_ENSG\\d+"), colnames(temp2), value = TRUE)
  
  if (length(gene_cols) > 1) {
    # If multiple columns for the gene, compute the row-wise average
    temp2[[gene]] <- rowMeans(temp2[gene_cols], na.rm = TRUE)
    # Remove the original ENSG columns after averaging
    temp2 <- temp2 %>% dplyr::select(-all_of(gene_cols))
  } else if (length(gene_cols) == 1) {
    # If only one column for the gene, just rename it by removing the ENSG part
    colnames(temp2)[colnames(temp2) == gene_cols] <- gene
  }
}

# View the result
print(temp2)

NA

temp.cor <- cor(temp2, method = "spearman")

p1 <- pheatmap(data.matrix(temp.cor), 
               scale = "none",
               cluster_rows = TRUE, 
               cluster_cols = TRUE,
               legend = TRUE,
               fontsize = 10)
p1

ggsave(filename = paste0(PLOT_loc, "/", Today, ".correlations.gene_list.pdf"), plot = p1, height = 15, width = 15)

ggsave(filename = paste0(PLOT_loc, "/", Today, ".correlations.gene_list.png"), plot = p1, height = 15, width = 15)


rm(temp, temp2, temp.cor, p1)

Gene target list

We are saving the final list of genes of interest

temp <- subset(expression_data_sel, select = c("gene_ensembl", "symbol"))

fwrite(temp,
       file = paste0(OUT_loc, "/", Today, ".target_list.qc.txt"),
       sep = " ", row.names = FALSE, col.names = TRUE,
       showProgress = TRUE)

Saving AERNA data subset

We will create a list of samples that should be included based on CEA, and having the proper informed consent (‘academic’). We will save the SummarizedExperiment as a RDS file for easy loading. The clinical data will also be saved as a separate txt-file.

To this end you need to lookup the EnsemblID (see above) for your target(s).

# targets = c("ENSG00000166033", # target
#             "ENSG00000169174") # positive control PCSK9
# targets = c("ENSG00000169174") # target/positive control PCSK9

targets = unlist(temp$gene_ensembl)

rm(temp)

Create a SummarizedExperiment

We make a SummarizedExperiment for the subset of RNAseq data. We make sure to only include the samples we need based on informed consent and we include only the requested variables.

Next, we are constructing the SummarizedExperiment.

AERNAS_meta_data = c("SampleType", "RNAseqTech", "RNAseqType", "RNAseqQC",
                     "StudyType", "StudyName", "Gender", "StudyBiobank", "SampleSize")

cat("\n* making a SummarizedExperiment ...\n")

* making a SummarizedExperiment ...

# this is all the data passing RNAseq quality control and UMI-corrected
# - after filtering on informed consent and artery type, the end sample size should be 1093
# - after filtering on 'no commercial business' based on informed consent, there are fewer samples: yyyy
cat("  > getting counts\n")

  > getting counts

bulkRNA_countsFilt <- as.data.frame(assay(AERNAScomboSE[targets,]))
dim(bulkRNA_countsFilt)

[1]   40 1092

head(bulkRNA_countsFilt)

cat("  > meta data and clinical data\n")

  > meta data and clinical data

# bulkRNA_meta_clin_COMMERCIAL <- subset(bulkRNA_meta_clin, select = c("study_number", basetable_vars))
bulkRNA_meta_clin_ACADEMIC <- subset(bulkRNA_meta_clin, select = c("study_number", 
                                                                   AERNAS_meta_data, # meta data of RNAseq experiment
                                                                   basetable_vars)) # selection of variables
dim(bulkRNA_meta_clin_ACADEMIC)

[1] 1092   68

head(bulkRNA_meta_clin_ACADEMIC)

cat("  > rowRanges\n")

  > rowRanges

# bulkRNA_meta_clin_COMMERCIAL <- subset(bulkRNA_meta_clin, select = c("study_number", basetable_vars))
bulkRNA_rowRanges <- rowRanges(AERNAScomboSE[targets])
bulkRNA_rowRanges

GRanges object with 40 ranges and 2 metadata columns:
                  seqnames              ranges strand |      feature_id      symbol
                     <Rle>           <IRanges>  <Rle> |     <character> <character>
  ENSG00000036530       14   99684298-99727301      + | ENSG00000036530     CYP46A1
  ENSG00000069869       15   55826922-55993660      - | ENSG00000069869       NEDD4
  ENSG00000080854       11 133896438-133956968      - | ENSG00000080854      IGSF9B
  ENSG00000087086       19   48965309-48966879      + | ENSG00000087086         FTL
  ENSG00000088827       20     3686970-3712600      - | ENSG00000088827     SIGLEC1
              ...      ...                 ...    ... .             ...         ...
  ENSG00000198734        1 169511951-169586588      - | ENSG00000198734          F5
  ENSG00000212907       MT         10470-10766      + | ENSG00000212907        ND4L
  ENSG00000213088        1 159203307-159206500      + | ENSG00000213088       ACKR1
  ENSG00000223417       11   89911111-89922245      - | ENSG00000223417    TRIM49D1
  ENSG00000228253       MT           8366-8572      + | ENSG00000228253        ATP8
  -------
  seqinfo: 331 sequences from an unspecified genome; no seqlengths

cat("  > construction of the SE\n")

  > construction of the SE

(AERNAScomboSE_target <- SummarizedExperiment(assays = list(counts = as.matrix(bulkRNA_countsFilt)),
                                              colData = bulkRNA_meta_clin_ACADEMIC,
                                              rowRanges = bulkRNA_rowRanges,
                                              metadata = "Athero-Express Biobank Study bulk RNA sequencing. Sample type: carotid plaques. Technology: CEL2-seq adapted for bulk RNA sequencing, thus 3'-focused. UMI-corrected. VST and size factor normalized. Target subset."))

class: RangedSummarizedExperiment 
dim: 40 1092 
metadata(1): ''
assays(1): counts
rownames(40): ENSG00000036530 ENSG00000069869 ... ENSG00000223417 ENSG00000228253
rowData names(2): feature_id symbol
colnames(1092): ae1 ae1026 ... ae986 ae992
colData names(68): study_number SampleType ... OverallPlaquePhenotype Plaque_Vulnerability_Index

cat("\n* removing intermediate files ...\n")

* removing intermediate files ...

rm(bulkRNA_meta_clin_ACADEMIC, bulkRNA_countsFilt, bulkRNA_rowRanges)

Saving data subset

Sample list and clinical data.

temp <- as.tibble(subset(colData(AERNAScomboSE_target), select = c("study_number", AERNAS_meta_data)))

fwrite(temp,
       file = paste0(OUT_loc, "/", Today, ".AERNAScomboSE_target.CEA.1092pts.samplelist.after_qc.IC_academic.csv"),
       sep = ",", row.names = FALSE, col.names = TRUE,
       showProgress = TRUE)
rm(temp)

temp <- as.tibble(subset(colData(AERNAScomboSE_target), select = c("study_number", basetable_vars)))

fwrite(temp,
       file = paste0(OUT_loc, "/", Today, ".AERNAScomboSE_target.CEA.1092pts.clinicaldata.after_qc.IC_academic.csv"),
       sep = ",", row.names = FALSE, col.names = TRUE,
       showProgress = TRUE)
rm(temp)

Assay data and gene information.

temp <- as_tibble(assay(AERNAScomboSE_target))

fwrite(temp,
       file = paste0(OUT_loc, "/", Today, ".AERNAScomboSE_target.CEA.1092pts.assay.after_qc.IC_academic.csv"),
       sep = ",", row.names = FALSE, col.names = TRUE,
       showProgress = TRUE)
rm(temp)

temp <- as_tibble(rowRanges(AERNAScomboSE_target))

fwrite(temp,
       file = paste0(OUT_loc, "/", Today, ".AERNAScomboSE_target.CEA.1092pts.genedata.after_qc.IC_academic.csv"),
       sep = ",", row.names = FALSE, col.names = TRUE,
       showProgress = TRUE)
rm(temp)

The subset as SummarizedExperiment() dataset in RDS-format.

# saveRDS(AERNASE, file = paste0(OUT_loc, "/", Today, ".AERNAScomboSE_target.CEA.608pts.SE.after_qc.IC_commercial.RDS"))
saveRDS(AERNAScomboSE_target, file = paste0(OUT_loc, "/", Today, ".AERNAScomboSE_target.CEA.1092pts.SE.after_qc.IC_academic.RDS"))

Plot selected data

Here we just do a sanity check and compare the expression for a favorite gene.

ggpubr::gghistogram(as.tibble(t(subset(assay(AERNAScomboSE_target), AERNAScomboSE_target@rowRanges$symbol == "HMOX1"))),
                    x = "ENSG00000100292", 
                    xlab = "HMOX1 (ENSG00000100292) expression\nlog-transformed, size corrected counts",
                    color = "white", fill =  uithof_color[20],
                    rug = F, add_density = F,
                    add = c("median"),
                    add.params = list(color = uithof_color[3]), 
                    ggtheme = theme_pubclean())

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.Warning: `geom_vline()`: Ignoring `mapping` because `xintercept` was provided.Warning: `geom_vline()`: Ignoring `data` because `xintercept` was provided.

ggpubr::gghistogram(as.tibble(t(subset(assay(AERNAScomboSE_target), AERNAScomboSE_target@rowRanges$symbol == "SORBS2"))),
                    x = "ENSG00000154556", 
                    xlab = "SORBS2 (ENSG00000154556) expression\nlog-transformed, size corrected counts",
                    color = "white", fill =  uithof_color[16],
                    rug = F, add_density = F,
                    add = c("median"),
                    add.params = list(color = uithof_color[3]), 
                    ggtheme = theme_pubclean())

Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.Warning: `geom_vline()`: Ignoring `mapping` because `xintercept` was provided.Warning: `geom_vline()`: Ignoring `data` because `xintercept` was provided.

Session information

---

Version:      v1.1.0
Last update:  2024-10-17
Written by:   Sander W. van der Laan (s.w.vanderlaan-2[at]umcutrecht.nl).
Description:  Script to load bulk RNA sequencing data, and perform gene expression analyses, and visualisations.
Minimum requirements: R version 3.5.2 (2018-12-20) -- 'Eggshell Igloo', macOS Mojave (10.14.2).

**MoSCoW To-Do List**
The things we Must, Should, Could, and Would have given the time we have.
_M_

_S_

_C_

_W_

**Changes log**
* v1.1.0 Added more detailed sub-analyses on the mean expression by gender. 
* v1.0.0 Inital version. Re-arranged to explore the subset of data. 

---

sessionInfo()

R version 4.4.1 (2024-06-14)
Platform: x86_64-apple-darwin20
Running under: macOS 15.1

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/New_York
tzcode source: internal

attached base packages:
 [1] stats4    grid      tools     stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] magrittr_2.0.3                          rmarkdown_2.28                          annotables_0.2.0                        EnhancedVolcano_1.22.0                 
 [5] ggrepel_0.9.6                           EnsDb.Hsapiens.v86_2.99.0               ensembldb_2.28.1                        AnnotationFilter_1.28.0                
 [9] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 mygene_1.40.0                           txdbmaker_1.0.1                         org.Hs.eg.db_3.19.1                    
[13] DESeq2_1.44.0                           SummarizedExperiment_1.34.0             MatrixGenerics_1.16.0                   matrixStats_1.4.1                      
[17] GenomicFeatures_1.56.0                  GenomicRanges_1.56.2                    AnnotationDbi_1.66.0                    Biobase_2.64.0                         
[21] GenomeInfoDb_1.40.1                     IRanges_2.38.1                          S4Vectors_0.42.1                        BiocGenerics_0.50.0                    
[25] Hmisc_5.1-3                             survminer_0.4.9                         survival_3.7-0                          GGally_2.2.1                           
[29] PerformanceAnalytics_2.0.4              xts_0.14.0                              zoo_1.8-12                              ggcorrplot_0.1.4.999                   
[33] corrr_0.4.4                             reshape2_1.4.4                          meta_7.0-0                              metadat_1.2-0                          
[37] qqman_0.1.9                             sjlabelled_1.2.0                        patchwork_1.3.0.9000                    labelled_2.13.0                        
[41] sjPlot_2.8.16                           UpSetR_1.4.0                            ggpubr_0.6.0.999                        forestplot_3.1.5                       
[45] abind_1.4-8                             checkmate_2.3.2                         pheatmap_1.0.12                         devtools_2.4.5                         
[49] usethis_3.0.0                           BlandAltmanLeh_0.3.1                    tableone_0.13.2                         openxlsx_4.2.7.1                       
[53] haven_2.5.4                             eeptools_1.2.5                          DT_0.33                                 knitr_1.48                             
[57] lubridate_1.9.3                         forcats_1.0.0                           stringr_1.5.1                           purrr_1.0.2                            
[61] tibble_3.2.1                            ggplot2_3.5.1                           tidyverse_2.0.0                         data.table_1.16.2                      
[65] naniar_1.1.0                            tidylog_1.1.0                           tidyr_1.3.1                             dplyr_1.1.4                            
[69] optparse_1.7.5                          readr_2.1.5                             pander_0.6.5                            R.utils_2.12.3                         
[73] R.oo_1.26.0                             R.methodsS3_1.8.2                       worcs_0.1.15                            credentials_2.0.2                      

loaded via a namespace (and not attached):
  [1] sjmisc_2.8.10            progress_1.2.3           urlchecker_1.0.1         nnet_7.3-19              Biostrings_2.72.1        vctrs_0.6.5             
  [7] rticles_0.27             digest_0.6.37            png_0.1-8                proxy_0.4-27             renv_1.0.11              MASS_7.3-61             
 [13] httpuv_1.6.15            withr_3.0.1              xfun_0.48                ellipsis_0.3.2           memoise_2.0.1            profvis_0.4.0           
 [19] ggsci_3.2.0              systemfonts_1.1.0        calibrate_1.7.7          ragg_1.3.3               Formula_1.2-5            sys_3.4.3               
 [25] prettyunits_1.2.0        datawizard_0.13.0        KEGGREST_1.44.1          promises_1.3.0           httr_1.4.7               rstatix_0.7.2           
 [31] restfulr_0.0.15          rstudioapi_0.16.0        UCSC.utils_1.0.0         miniUI_0.1.1.1           generics_0.1.3           base64enc_0.1-3         
 [37] curl_5.2.3               mitools_2.4              zlibbioc_1.50.0          ggeffects_1.7.2          GenomeInfoDbData_1.2.12  quadprog_1.5-8          
 [43] SparseArray_1.4.8        xtable_1.8-4             evaluate_1.0.1           S4Arrays_1.4.1           BiocFileCache_2.12.0     hms_1.1.3               
 [49] colorspace_2.1-1         filelock_1.0.3           getopt_1.20.4            lmtest_0.9-40            later_1.3.2              lattice_0.22-6          
 [55] XML_3.99-0.17            survey_4.4-2             class_7.3-22             pillar_1.9.0             nlme_3.1-166             performance_0.12.3      
 [61] compiler_4.4.1           stringi_1.8.4            minqa_1.2.8              GenomicAlignments_1.40.0 plyr_1.8.9               crayon_1.5.3            
 [67] BiocIO_1.14.0            chron_2.3-61             locfit_1.5-9.10          bit_4.5.0                mathjaxr_1.6-0           textshaping_0.4.0       
 [73] codetools_0.2-20         clisymbols_1.2.0         openssl_2.2.2            crosstalk_1.2.1          bslib_0.8.0              e1071_1.7-16            
 [79] mime_0.12                splines_4.4.1            Rcpp_1.0.13              dbplyr_2.5.0             blob_1.2.4               utf8_1.2.4              
 [85] lme4_1.1-35.5            fs_1.6.4                 prereg_0.6.0             pkgbuild_1.4.4           gh_1.4.1                 ggsignif_0.6.4          
 [91] sqldf_0.4-11             Matrix_1.7-0             tzdb_0.4.0               visdat_0.6.0             pkgconfig_2.0.3          cachem_1.1.0            
 [97] RSQLite_2.3.7            DBI_1.2.3                numDeriv_2016.8-1.1      fastmap_1.2.0            scales_1.3.0             Rsamtools_2.20.0        
[103] broom_1.0.7              sass_0.4.9               coda_0.19-4.1            ggstats_0.7.0            insight_0.20.5           carData_3.0-5           
[109] rpart_4.1.23             farver_2.1.2             gsubfn_0.7               yaml_2.3.10              foreign_0.8-87           sjstats_0.19.0          
[115] rtracklayer_1.64.0       cli_3.6.3                lifecycle_1.0.4          askpass_1.2.1            sessioninfo_1.2.2        backports_1.5.0         
[121] BiocParallel_1.38.0      timechange_0.3.0         gtable_0.3.5             rjson_0.2.23             metafor_4.6-0            parallel_4.4.1          
[127] jsonlite_1.8.9           bitops_1.0-9             bit64_4.5.2              proto_1.0.0              zip_2.3.1                ranger_0.16.0           
[133] jquerylib_0.1.4          survMisc_0.5.6           lazyeval_0.2.2           shiny_1.9.1              htmltools_0.5.8.1        KMsurv_0.1-5            
[139] rappdirs_0.3.3           tinytex_0.53             glue_1.8.0               httr2_1.0.5              XVector_0.44.0           RCurl_1.98-1.16         
[145] gridExtra_2.3            boot_1.3-31              R6_2.5.1                 arm_1.14-4               km.ci_0.5-6              vcd_1.4-13              
[151] CompQuadForm_1.4.3       labeling_0.4.3           cluster_2.1.6            pkgload_1.4.0            nloptr_2.1.1             DelayedArray_0.30.1     
[157] tidyselect_1.2.1         ProtGenerics_1.36.0      htmlTable_2.4.3          xml2_1.3.6               car_3.1-3                munsell_0.5.1           
[163] htmlwidgets_1.6.4        RColorBrewer_1.1-3       biomaRt_2.60.1           rlang_1.1.4              gert_2.1.4               remotes_2.5.0           
[169] fansi_1.0.6             

Saving environment

save.image(paste0(PROJECT_loc, "/",Today,".",PROJECTNAME,".bulkRNAseq.exploration.RData"))

© 1979-2024 Sander W. van der Laan | s.w.vanderlaan[at]gmail[dot]com | vanderlaanand.science.